acc1_cancer_cells = readRDS("./Data/acc1_cancer_cells_15KnCount_V3.RDS")
all_acc_cancer_cells = readRDS("./Data/acc_cancer_cells_V4.RDS")
acc_all_cells = readRDS("./Data/acc_tpm_nCount_mito_no146_15k_with_ACC1_.RDS")
acc_cancerCells_noACC1 = readRDS("./Data/acc_cancer_no146_primaryonly15k_cancercells.rds")
luminal_pathways = c("CHARAFE_BREAST_CANCER_LUMINAL_VS_BASAL_DN","CHARAFE_BREAST_CANCER_LUMINAL_VS_BASAL_UP","CHARAFE_BREAST_CANCER_LUMINAL_VS_MESENCHYMAL_DN","CHARAFE_BREAST_CANCER_LUMINAL_VS_MESENCHYMAL_UP","HUPER_BREAST_BASAL_VS_LUMINAL_DN","LIM_MAMMARY_LUMINAL_PROGENITOR_UP","SMID_BREAST_CANCER_LUMINAL_B_UP" )
# add luminal pathwaysx
luminal_gs = msigdbr(species = "Homo sapiens") %>%as.data.frame() %>% dplyr::filter(gs_name %in% luminal_pathways)%>% dplyr::distinct(gs_name, gene_symbol) %>% as.data.frame()
source_from_github(repositoy = "DEG_functions",version = "0.2.24")
source_from_github(repositoy = "HMSC_functions",version = "0.1.13",script_name = "functions.R")
source_from_github(repositoy = "cNMF_functions",version = "0.3.72",script_name = "cnmf_function_Harmony.R")
DimPlot(object = acc1_cancer_cells,pt.size = 2)
deg = FindAllMarkers(object = acc1_cancer_cells,features = VariableFeatures(acc1_cancer_cells),densify = T)
Calculating cluster 0
| | 0 % ~calculating
|+ | 1 % ~04s
|+ | 2 % ~04s
|++ | 3 % ~03s
|++ | 4 % ~03s
|+++ | 5 % ~03s
|+++ | 6 % ~03s
|++++ | 7 % ~03s
|++++ | 8 % ~03s
|+++++ | 9 % ~03s
|+++++ | 10% ~03s
|++++++ | 11% ~03s
|++++++ | 12% ~03s
|+++++++ | 13% ~03s
|+++++++ | 14% ~03s
|++++++++ | 15% ~03s
|++++++++ | 16% ~03s
|+++++++++ | 17% ~03s
|+++++++++ | 18% ~03s
|++++++++++ | 19% ~03s
|++++++++++ | 20% ~03s
|+++++++++++ | 21% ~03s
|+++++++++++ | 22% ~03s
|++++++++++++ | 23% ~03s
|++++++++++++ | 24% ~03s
|+++++++++++++ | 25% ~02s
|+++++++++++++ | 26% ~02s
|++++++++++++++ | 27% ~02s
|++++++++++++++ | 28% ~02s
|+++++++++++++++ | 29% ~02s
|+++++++++++++++ | 30% ~02s
|++++++++++++++++ | 31% ~02s
|++++++++++++++++ | 32% ~03s
|+++++++++++++++++ | 33% ~03s
|+++++++++++++++++ | 34% ~02s
|++++++++++++++++++ | 35% ~02s
|++++++++++++++++++ | 36% ~02s
|+++++++++++++++++++ | 37% ~02s
|+++++++++++++++++++ | 38% ~02s
|++++++++++++++++++++ | 39% ~02s
|++++++++++++++++++++ | 40% ~02s
|+++++++++++++++++++++ | 41% ~02s
|+++++++++++++++++++++ | 42% ~02s
|++++++++++++++++++++++ | 43% ~02s
|++++++++++++++++++++++ | 44% ~02s
|+++++++++++++++++++++++ | 45% ~02s
|+++++++++++++++++++++++ | 46% ~02s
|++++++++++++++++++++++++ | 47% ~02s
|++++++++++++++++++++++++ | 48% ~02s
|+++++++++++++++++++++++++ | 49% ~02s
|+++++++++++++++++++++++++ | 50% ~02s
|++++++++++++++++++++++++++ | 51% ~02s
|++++++++++++++++++++++++++ | 52% ~02s
|+++++++++++++++++++++++++++ | 53% ~02s
|+++++++++++++++++++++++++++ | 54% ~02s
|++++++++++++++++++++++++++++ | 55% ~02s
|++++++++++++++++++++++++++++ | 56% ~02s
|+++++++++++++++++++++++++++++ | 57% ~02s
|+++++++++++++++++++++++++++++ | 58% ~02s
|++++++++++++++++++++++++++++++ | 59% ~02s
|++++++++++++++++++++++++++++++ | 60% ~01s
|+++++++++++++++++++++++++++++++ | 61% ~01s
|+++++++++++++++++++++++++++++++ | 62% ~01s
|++++++++++++++++++++++++++++++++ | 63% ~01s
|++++++++++++++++++++++++++++++++ | 64% ~01s
|+++++++++++++++++++++++++++++++++ | 65% ~01s
|+++++++++++++++++++++++++++++++++ | 66% ~01s
|++++++++++++++++++++++++++++++++++ | 67% ~01s
|++++++++++++++++++++++++++++++++++ | 68% ~01s
|+++++++++++++++++++++++++++++++++++ | 69% ~01s
|+++++++++++++++++++++++++++++++++++ | 70% ~01s
|++++++++++++++++++++++++++++++++++++ | 71% ~01s
|++++++++++++++++++++++++++++++++++++ | 72% ~01s
|+++++++++++++++++++++++++++++++++++++ | 73% ~01s
|+++++++++++++++++++++++++++++++++++++ | 74% ~01s
|++++++++++++++++++++++++++++++++++++++ | 75% ~01s
|++++++++++++++++++++++++++++++++++++++ | 76% ~01s
|+++++++++++++++++++++++++++++++++++++++ | 77% ~01s
|+++++++++++++++++++++++++++++++++++++++ | 78% ~01s
|++++++++++++++++++++++++++++++++++++++++ | 79% ~01s
|++++++++++++++++++++++++++++++++++++++++ | 80% ~01s
|+++++++++++++++++++++++++++++++++++++++++ | 81% ~01s
|+++++++++++++++++++++++++++++++++++++++++ | 82% ~01s
|++++++++++++++++++++++++++++++++++++++++++ | 83% ~01s
|++++++++++++++++++++++++++++++++++++++++++ | 84% ~01s
|+++++++++++++++++++++++++++++++++++++++++++ | 85% ~01s
|+++++++++++++++++++++++++++++++++++++++++++ | 86% ~00s
|++++++++++++++++++++++++++++++++++++++++++++ | 87% ~00s
|++++++++++++++++++++++++++++++++++++++++++++ | 88% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++ | 89% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++ | 90% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++ | 91% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++ | 92% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++ | 93% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++ | 95% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++ | 96% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=04s
Calculating cluster 1
| | 0 % ~calculating
|+ | 1 % ~04s
|+ | 2 % ~04s
|++ | 3 % ~05s
|++ | 4 % ~04s
|+++ | 5 % ~04s
|+++ | 6 % ~04s
|++++ | 7 % ~04s
|++++ | 8 % ~04s
|+++++ | 9 % ~03s
|+++++ | 10% ~03s
|++++++ | 11% ~03s
|++++++ | 12% ~03s
|+++++++ | 13% ~03s
|+++++++ | 14% ~03s
|++++++++ | 15% ~03s
|++++++++ | 16% ~03s
|+++++++++ | 17% ~03s
|+++++++++ | 18% ~03s
|++++++++++ | 19% ~03s
|++++++++++ | 20% ~03s
|+++++++++++ | 21% ~03s
|+++++++++++ | 22% ~03s
|++++++++++++ | 23% ~03s
|++++++++++++ | 24% ~03s
|+++++++++++++ | 25% ~03s
|+++++++++++++ | 26% ~03s
|++++++++++++++ | 27% ~02s
|++++++++++++++ | 28% ~02s
|+++++++++++++++ | 29% ~02s
|+++++++++++++++ | 30% ~02s
|++++++++++++++++ | 31% ~02s
|++++++++++++++++ | 32% ~02s
|+++++++++++++++++ | 33% ~02s
|+++++++++++++++++ | 34% ~02s
|++++++++++++++++++ | 35% ~02s
|++++++++++++++++++ | 36% ~02s
|+++++++++++++++++++ | 37% ~02s
|+++++++++++++++++++ | 38% ~02s
|++++++++++++++++++++ | 39% ~02s
|++++++++++++++++++++ | 40% ~02s
|+++++++++++++++++++++ | 41% ~02s
|+++++++++++++++++++++ | 42% ~02s
|++++++++++++++++++++++ | 43% ~02s
|++++++++++++++++++++++ | 44% ~02s
|+++++++++++++++++++++++ | 45% ~02s
|+++++++++++++++++++++++ | 46% ~02s
|++++++++++++++++++++++++ | 47% ~02s
|++++++++++++++++++++++++ | 48% ~02s
|+++++++++++++++++++++++++ | 49% ~02s
|+++++++++++++++++++++++++ | 50% ~02s
|++++++++++++++++++++++++++ | 51% ~02s
|++++++++++++++++++++++++++ | 52% ~02s
|+++++++++++++++++++++++++++ | 53% ~02s
|+++++++++++++++++++++++++++ | 54% ~02s
|++++++++++++++++++++++++++++ | 55% ~02s
|++++++++++++++++++++++++++++ | 56% ~01s
|+++++++++++++++++++++++++++++ | 57% ~01s
|+++++++++++++++++++++++++++++ | 58% ~02s
|++++++++++++++++++++++++++++++ | 59% ~02s
|++++++++++++++++++++++++++++++ | 60% ~01s
|+++++++++++++++++++++++++++++++ | 61% ~01s
|+++++++++++++++++++++++++++++++ | 62% ~01s
|++++++++++++++++++++++++++++++++ | 63% ~01s
|++++++++++++++++++++++++++++++++ | 64% ~01s
|+++++++++++++++++++++++++++++++++ | 65% ~01s
|+++++++++++++++++++++++++++++++++ | 66% ~01s
|++++++++++++++++++++++++++++++++++ | 67% ~01s
|++++++++++++++++++++++++++++++++++ | 68% ~01s
|+++++++++++++++++++++++++++++++++++ | 69% ~01s
|+++++++++++++++++++++++++++++++++++ | 70% ~01s
|++++++++++++++++++++++++++++++++++++ | 71% ~01s
|++++++++++++++++++++++++++++++++++++ | 72% ~01s
|+++++++++++++++++++++++++++++++++++++ | 73% ~01s
|+++++++++++++++++++++++++++++++++++++ | 74% ~01s
|++++++++++++++++++++++++++++++++++++++ | 75% ~01s
|++++++++++++++++++++++++++++++++++++++ | 76% ~01s
|+++++++++++++++++++++++++++++++++++++++ | 77% ~01s
|+++++++++++++++++++++++++++++++++++++++ | 78% ~01s
|++++++++++++++++++++++++++++++++++++++++ | 79% ~01s
|++++++++++++++++++++++++++++++++++++++++ | 80% ~01s
|+++++++++++++++++++++++++++++++++++++++++ | 81% ~01s
|+++++++++++++++++++++++++++++++++++++++++ | 82% ~01s
|++++++++++++++++++++++++++++++++++++++++++ | 83% ~01s
|++++++++++++++++++++++++++++++++++++++++++ | 84% ~01s
|+++++++++++++++++++++++++++++++++++++++++++ | 85% ~01s
|+++++++++++++++++++++++++++++++++++++++++++ | 86% ~01s
|++++++++++++++++++++++++++++++++++++++++++++ | 87% ~00s
|++++++++++++++++++++++++++++++++++++++++++++ | 88% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++ | 89% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++ | 90% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++ | 91% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++ | 92% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++ | 93% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++ | 95% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++ | 96% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=04s
Calculating cluster 2
| | 0 % ~calculating
|+ | 1 % ~04s
|++ | 2 % ~04s
|++ | 3 % ~04s
|+++ | 4 % ~04s
|+++ | 5 % ~03s
|++++ | 6 % ~03s
|++++ | 7 % ~03s
|+++++ | 8 % ~03s
|+++++ | 9 % ~03s
|++++++ | 10% ~03s
|++++++ | 11% ~03s
|+++++++ | 12% ~03s
|+++++++ | 13% ~03s
|++++++++ | 14% ~03s
|++++++++ | 15% ~03s
|+++++++++ | 16% ~03s
|+++++++++ | 17% ~03s
|++++++++++ | 18% ~03s
|++++++++++ | 19% ~03s
|+++++++++++ | 20% ~03s
|+++++++++++ | 21% ~03s
|++++++++++++ | 22% ~03s
|++++++++++++ | 23% ~03s
|+++++++++++++ | 24% ~03s
|+++++++++++++ | 25% ~03s
|++++++++++++++ | 26% ~03s
|++++++++++++++ | 27% ~03s
|+++++++++++++++ | 28% ~03s
|+++++++++++++++ | 29% ~03s
|++++++++++++++++ | 30% ~03s
|++++++++++++++++ | 31% ~02s
|+++++++++++++++++ | 32% ~02s
|+++++++++++++++++ | 33% ~02s
|++++++++++++++++++ | 34% ~02s
|++++++++++++++++++ | 35% ~02s
|+++++++++++++++++++ | 36% ~02s
|+++++++++++++++++++ | 37% ~02s
|++++++++++++++++++++ | 38% ~02s
|++++++++++++++++++++ | 39% ~02s
|+++++++++++++++++++++ | 40% ~02s
|+++++++++++++++++++++ | 41% ~02s
|++++++++++++++++++++++ | 42% ~02s
|++++++++++++++++++++++ | 43% ~02s
|+++++++++++++++++++++++ | 44% ~02s
|+++++++++++++++++++++++ | 45% ~02s
|++++++++++++++++++++++++ | 46% ~02s
|++++++++++++++++++++++++ | 47% ~02s
|+++++++++++++++++++++++++ | 48% ~02s
|+++++++++++++++++++++++++ | 49% ~02s
|++++++++++++++++++++++++++ | 51% ~02s
|++++++++++++++++++++++++++ | 52% ~02s
|+++++++++++++++++++++++++++ | 53% ~02s
|+++++++++++++++++++++++++++ | 54% ~02s
|++++++++++++++++++++++++++++ | 55% ~02s
|++++++++++++++++++++++++++++ | 56% ~02s
|+++++++++++++++++++++++++++++ | 57% ~01s
|+++++++++++++++++++++++++++++ | 58% ~01s
|++++++++++++++++++++++++++++++ | 59% ~01s
|++++++++++++++++++++++++++++++ | 60% ~01s
|+++++++++++++++++++++++++++++++ | 61% ~01s
|+++++++++++++++++++++++++++++++ | 62% ~01s
|++++++++++++++++++++++++++++++++ | 63% ~01s
|++++++++++++++++++++++++++++++++ | 64% ~01s
|+++++++++++++++++++++++++++++++++ | 65% ~01s
|+++++++++++++++++++++++++++++++++ | 66% ~01s
|++++++++++++++++++++++++++++++++++ | 67% ~01s
|++++++++++++++++++++++++++++++++++ | 68% ~01s
|+++++++++++++++++++++++++++++++++++ | 69% ~01s
|+++++++++++++++++++++++++++++++++++ | 70% ~01s
|++++++++++++++++++++++++++++++++++++ | 71% ~01s
|++++++++++++++++++++++++++++++++++++ | 72% ~01s
|+++++++++++++++++++++++++++++++++++++ | 73% ~01s
|+++++++++++++++++++++++++++++++++++++ | 74% ~01s
|++++++++++++++++++++++++++++++++++++++ | 75% ~01s
|++++++++++++++++++++++++++++++++++++++ | 76% ~01s
|+++++++++++++++++++++++++++++++++++++++ | 77% ~01s
|+++++++++++++++++++++++++++++++++++++++ | 78% ~01s
|++++++++++++++++++++++++++++++++++++++++ | 79% ~01s
|++++++++++++++++++++++++++++++++++++++++ | 80% ~01s
|+++++++++++++++++++++++++++++++++++++++++ | 81% ~01s
|+++++++++++++++++++++++++++++++++++++++++ | 82% ~01s
|++++++++++++++++++++++++++++++++++++++++++ | 83% ~01s
|++++++++++++++++++++++++++++++++++++++++++ | 84% ~01s
|+++++++++++++++++++++++++++++++++++++++++++ | 85% ~01s
|+++++++++++++++++++++++++++++++++++++++++++ | 86% ~01s
|++++++++++++++++++++++++++++++++++++++++++++ | 87% ~00s
|++++++++++++++++++++++++++++++++++++++++++++ | 88% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++ | 89% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++ | 90% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++ | 91% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++ | 92% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++ | 93% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++ | 95% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++ | 96% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=04s
Calculating cluster 3
| | 0 % ~calculating
|+ | 1 % ~04s
|+ | 2 % ~04s
|++ | 3 % ~03s
|++ | 4 % ~03s
|+++ | 5 % ~03s
|+++ | 6 % ~03s
|++++ | 7 % ~03s
|++++ | 8 % ~03s
|+++++ | 9 % ~03s
|+++++ | 10% ~03s
|++++++ | 11% ~03s
|++++++ | 12% ~03s
|+++++++ | 13% ~03s
|+++++++ | 14% ~03s
|++++++++ | 15% ~03s
|++++++++ | 16% ~03s
|+++++++++ | 17% ~03s
|+++++++++ | 18% ~03s
|++++++++++ | 19% ~03s
|++++++++++ | 20% ~03s
|+++++++++++ | 21% ~03s
|+++++++++++ | 22% ~03s
|++++++++++++ | 23% ~03s
|++++++++++++ | 24% ~03s
|+++++++++++++ | 25% ~03s
|+++++++++++++ | 26% ~03s
|++++++++++++++ | 27% ~03s
|++++++++++++++ | 28% ~02s
|+++++++++++++++ | 29% ~02s
|+++++++++++++++ | 30% ~02s
|++++++++++++++++ | 31% ~02s
|++++++++++++++++ | 32% ~02s
|+++++++++++++++++ | 33% ~02s
|+++++++++++++++++ | 34% ~02s
|++++++++++++++++++ | 35% ~02s
|++++++++++++++++++ | 36% ~02s
|+++++++++++++++++++ | 37% ~02s
|+++++++++++++++++++ | 38% ~02s
|++++++++++++++++++++ | 39% ~02s
|++++++++++++++++++++ | 40% ~02s
|+++++++++++++++++++++ | 41% ~02s
|+++++++++++++++++++++ | 42% ~02s
|++++++++++++++++++++++ | 43% ~02s
|++++++++++++++++++++++ | 44% ~02s
|+++++++++++++++++++++++ | 45% ~02s
|+++++++++++++++++++++++ | 46% ~02s
|++++++++++++++++++++++++ | 47% ~02s
|++++++++++++++++++++++++ | 48% ~02s
|+++++++++++++++++++++++++ | 49% ~02s
|+++++++++++++++++++++++++ | 50% ~02s
|++++++++++++++++++++++++++ | 51% ~02s
|++++++++++++++++++++++++++ | 52% ~02s
|+++++++++++++++++++++++++++ | 53% ~02s
|+++++++++++++++++++++++++++ | 54% ~02s
|++++++++++++++++++++++++++++ | 55% ~02s
|++++++++++++++++++++++++++++ | 56% ~02s
|+++++++++++++++++++++++++++++ | 57% ~02s
|+++++++++++++++++++++++++++++ | 58% ~01s
|++++++++++++++++++++++++++++++ | 59% ~01s
|++++++++++++++++++++++++++++++ | 60% ~01s
|+++++++++++++++++++++++++++++++ | 61% ~01s
|+++++++++++++++++++++++++++++++ | 62% ~01s
|++++++++++++++++++++++++++++++++ | 63% ~01s
|++++++++++++++++++++++++++++++++ | 64% ~01s
|+++++++++++++++++++++++++++++++++ | 65% ~01s
|+++++++++++++++++++++++++++++++++ | 66% ~01s
|++++++++++++++++++++++++++++++++++ | 67% ~01s
|++++++++++++++++++++++++++++++++++ | 68% ~01s
|+++++++++++++++++++++++++++++++++++ | 69% ~01s
|+++++++++++++++++++++++++++++++++++ | 70% ~01s
|++++++++++++++++++++++++++++++++++++ | 71% ~01s
|++++++++++++++++++++++++++++++++++++ | 72% ~01s
|+++++++++++++++++++++++++++++++++++++ | 73% ~01s
|+++++++++++++++++++++++++++++++++++++ | 74% ~01s
|++++++++++++++++++++++++++++++++++++++ | 75% ~01s
|++++++++++++++++++++++++++++++++++++++ | 76% ~01s
|+++++++++++++++++++++++++++++++++++++++ | 77% ~01s
|+++++++++++++++++++++++++++++++++++++++ | 78% ~01s
|++++++++++++++++++++++++++++++++++++++++ | 79% ~01s
|++++++++++++++++++++++++++++++++++++++++ | 80% ~01s
|+++++++++++++++++++++++++++++++++++++++++ | 81% ~01s
|+++++++++++++++++++++++++++++++++++++++++ | 82% ~01s
|++++++++++++++++++++++++++++++++++++++++++ | 83% ~01s
|++++++++++++++++++++++++++++++++++++++++++ | 84% ~01s
|+++++++++++++++++++++++++++++++++++++++++++ | 85% ~01s
|+++++++++++++++++++++++++++++++++++++++++++ | 86% ~00s
|++++++++++++++++++++++++++++++++++++++++++++ | 87% ~00s
|++++++++++++++++++++++++++++++++++++++++++++ | 88% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++ | 89% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++ | 90% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++ | 91% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++ | 92% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++ | 93% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++ | 95% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++ | 96% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s
for (cluster in unique(deg$cluster)) {
print(deg[deg$cluster == cluster,])
}
for (cluster in unique(deg$cluster)) {
deg_of_cluster = deg[deg$cluster == cluster,]
enrichment_analysis(deg_of_cluster,background =VariableFeatures(acc1_cancer_cells),fdr_Cutoff = 0.01,ident.1 = paste("Cluster",cluster),ident.2 = paste("Cluster",cluster),show_by = 1)
}
NA
NA
FeaturePlot(object = acc1_cancer_cells,features = c("TP63","ACTA2","IL12B","CNN1"))
Warning in FeaturePlot(object = acc1_cancer_cells, features = c("TP63", :
All cells have the same value (0) of IL12B.
all_acc_cancer_cells = SetIdent(all_acc_cancer_cells, value ="patient.ident")
acc_deg <- FindMarkers(all_acc_cancer_cells, ident.1 = "HMSC",logfc.threshold = 1.5,features = VariableFeatures(all_acc_cancer_cells))
enrichment_analysis(acc_deg,background = VariableFeatures(all_acc_cancer_cells),fdr_Cutoff = 0.01,ident.1 = "HMSC",ident.2 = "ACC",show_by = 1)
hallmark_name = "GO_MITOTIC_CELL_CYCLE"
genesets =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
var_features=acc_cancerCells_noACC1@assays$RNA@var.features
geneIds= genesets[[hallmark_name]]@geneIds
score <- apply(acc_cancerCells_noACC1@assays$RNA@scale.data[intersect(geneIds,var_features),],2,mean)
acc_cancerCells_noACC1=AddMetaData(acc_cancerCells_noACC1,score,hallmark_name)
hallmark_name = "GO_MITOTIC_CELL_CYCLE"
genesets =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
var_features=acc1_cancer_cells@assays$RNA@var.features
geneIds= genesets[[hallmark_name]]@geneIds
score <- apply(acc1_cancer_cells@assays$RNA@scale.data[intersect(geneIds,var_features),],2,mean)
acc1_cancer_cells=AddMetaData(acc1_cancer_cells,score,hallmark_name)
# hallmark_name = "GO_MITOTIC_CELL_CYCLE"
# genesets =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
# var_features=acc1_cancer_cells@assays$RNA@var.features
# geneIds= genesets[[hallmark_name]]@geneIds
# score <- apply(acc1_cancer_cells@assays$RNA@scale.data[intersect(geneIds,var_features),],2,mean)
# acc1_cancer_cells=AddMetaData(acc1_cancer_cells,score,hallmark_name)
library(highcharter)
options(highcharter.theme = hc_theme_smpl(tooltip = list(valueDecimals = 2)))
cc_scores = FetchData(object = acc1_cancer_cells,vars = "GO_MITOTIC_CELL_CYCLE")
hchart(
density(cc_scores$GO_MITOTIC_CELL_CYCLE),
type = "area", name = "GO_MITOTIC_CELL_CYCLE"
)
NA
cc_scores = FetchData(object = acc_cancerCells_noACC1,vars = "GO_MITOTIC_CELL_CYCLE")
hchart(
density(cc_scores$GO_MITOTIC_CELL_CYCLE),
type = "area", name = "GO_MITOTIC_CELL_CYCLE"
)
NA
acc_cc_scores = FetchData(object = acc_cancerCells_noACC1,vars = "GO_MITOTIC_CELL_CYCLE")
hmsc_cc_scores = FetchData(object = acc1_cancer_cells,vars = "GO_MITOTIC_CELL_CYCLE")
ggplot() +
geom_density(aes(GO_MITOTIC_CELL_CYCLE, fill = "ACC"), alpha = .2, data = acc_cc_scores) +
geom_density(aes(GO_MITOTIC_CELL_CYCLE, fill = "HMSC"), alpha = .2, data = hmsc_cc_scores) +
scale_fill_manual(name = "Dataset", values = c(ACC = "red", HMSC = "green"))+ geom_vline(aes(xintercept=0.15),
color="blue", linetype="dashed", size=1)
#filter:
all_acc_cancer_cells_ccFiltered=all_acc_cancer_cells[,all_acc_cancer_cells@meta.data[[hallmark_name]]< 0.15]
min_threshold = min(all_acc_cancer_cells$GO_MITOTIC_CELL_CYCLE)
max_threshold = max(all_acc_cancer_cells$GO_MITOTIC_CELL_CYCLE)
library(viridis)
Loading required package: viridisLite
FeaturePlot(object = all_acc_cancer_cells,features = hallmark_name) + ggtitle("Before cc filtering") & scale_color_gradientn(colours = plasma(n = 10, direction = -1), limits = c(min_threshold, max_threshold))
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Warning in grSoftVersion() :
unable to load shared object '/usr/local/lib/R/modules//R_X11.so':
libXt.so.6: cannot open shared object file: No such file or directory
FeaturePlot(object = all_acc_cancer_cells_ccFiltered,features = hallmark_name) + ggtitle("After cc filtering") & scale_color_gradientn(colours = plasma(n = 10, direction = -1), limits = c(min_threshold, max_threshold))
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
1
patient.ident = all_acc_cancer_cells_ccFiltered$patient.ident %>% as.data.frame()
patient.ident[,1] = as.character(patient.ident[,1])
patient.ident[patient.ident[,1] == "ACC1",] = "HMSC"
patient.ident[,1] = as.factor(patient.ident[,1])
all_acc_cancer_cells_ccFiltered = AddMetaData(object = all_acc_cancer_cells_ccFiltered,metadata = patient.ident,col.name = "patient.ident")
all_acc_cancer_cells_ccFiltered = SetIdent(all_acc_cancer_cells_ccFiltered, value ="patient.ident")
acc_deg <- FindMarkers(all_acc_cancer_cells_ccFiltered, ident.1 = "HMSC",logfc.threshold = 1.5,features = VariableFeatures(all_acc_cancer_cells_ccFiltered))
enrichment_analysis(acc_deg,background = VariableFeatures(all_acc_cancer_cells_ccFiltered),fdr_Cutoff = 0.01,ident.1 = "HMSC",ident.2 = "ACC",show_by = 1)
all_acc_cancer_cells = SetIdent(object = all_acc_cancer_cells,value = "patient.ident") #active snn graph
FeaturePlot(object = all_acc_cancer_cells,features = "MYB",label = T)
all_acc_cancer_cells = SetIdent(object = all_acc_cancer_cells,value = "patient.ident") #active snn graph
FeaturePlot(object = all_acc_cancer_cells,features = "MYB",label = T)
new.cluster.ids <- c("cancer", #0
"cancer", #1
"CAF", #2
"cancer", #3
"Endothelial", #4
"cancer", #5
"cancer", #6
"CAF", #7
"CAF", #8
"CAF", #9
"cancer", #10
"CAF", #11
"cancer", #12
"cancer", #13
"cancer", #14
"cancer", #15
"cancer", #16
"WBC", #17
"CAF" #18
)
#rename idents:
acc_all_cells = SetIdent(object = acc_all_cells,value = "RNA_snn_res.1") #active snn graph
names(new.cluster.ids) <- levels(acc_all_cells) #add snn graph levels to new.cluster.ids
acc_all_cells@meta.data[["seurat_clusters"]] = acc_all_cells@meta.data[["RNA_snn_res.1"]]
acc_all_cells = SetIdent(object = acc_all_cells,value = "seurat_clusters")
acc_all_cells <- RenameIdents(acc_all_cells, new.cluster.ids)
# divide "cancer" into patients:
cell_types = acc_all_cells@active.ident %>% as.data.frame()
cell_types[,1]<- as.character(cell_types[,1])
cell_types = cbind(cell_types,acc_all_cells$patient.ident) %>% setnames(old = names(.),
new = c('cell_type','patient'))
cell_types[cell_types$cell_type == "cancer",] = cell_types[cell_types$cell_type == "cancer",2]
# hmsc_rows = (startsWith(x = rownames(cell_types),prefix = "ACC.plate2") | startsWith(x = rownames(cell_types),prefix = "ACC1.")) & cell_types[,1] == "cancer"
# acc_rows = !(startsWith(x = rownames(cell_types),prefix = "ACC.plate2") | startsWith(x = rownames(cell_types),prefix = "ACC1.")) & cell_types[,1] == "cancer"
# cell_types[,1][hmsc_rows] = "HMSC"
# cell_types[,1][acc_rows] = "ACC"
#add to metadata:
cell_types[,2] = NULL
cell_types[cell_types$cell_type == "ACC1",] = "HMSC"
acc_all_cells = AddMetaData(object =acc_all_cells ,metadata = cell_types,col.name = "cell.type")
library(infercnv)
library(tidyverse)
acc_annotation = as.data.frame(acc_all_cells@meta.data[,"cell.type",drop = F])
acc_annotation = acc_annotation %>% rownames_to_column("orig.ident")
acc_annotation = acc_annotation %>% mutate(orig.ident = gsub(x = acc_annotation$orig.ident,pattern = "\\.", replacement = "-") %>%
gsub(pattern = "_", replacement = "-"))
write.table(acc_annotation, "./Data/inferCNV/acc_annotation.txt", append = FALSE,
sep = "\t", dec = ".",row.names = FALSE, col.names = F)
infercnv_obj = CreateInfercnvObject(raw_counts_matrix="./Data/inferCNV/all.4icnv.txt",
annotations_file="./Data/inferCNV/acc_annotation.txt",
delim="\t",gene_order_file="./Data/inferCNV/gencode_v19_gene_pos.txt"
,ref_group_names=c("CAF", "Endothelial", "WBC")) #groups of normal cells
infercnv_obj_default = infercnv::run(infercnv_obj, cutoff=1, out_dir='./Data/inferCNV/infercnv_output',
cluster_by_groups=T, plot_steps=FALSE,
denoise=TRUE, HMM=FALSE, no_prelim_plot=TRUE,
png_res=300)
plot_cnv(infercnv_obj_default, output_format = "png", write_expr_matrix = FALSE,out_dir = "./Data/inferCNV/",png_res =800,obs_title = "Malignant cells",ref_title = "Normal cells")
cluster.info=FetchData(acc_all_cells,c("ident","orig.ident","UMAP_1","UMAP_2","nCount_RNA","nFeature_RNA","percent.mt","patient.ident","seurat_clusters"))
cluster.info$cell=rownames(cluster.info)
library(limma)
smoothed=apply(infercnv_obj_default@expr.data,2,tricubeMovingAverage, span=0.01)
cnsig=sqrt(apply((smoothed-1)^2,2,mean))
umap=cluster.info
names(cnsig)=umap$cell
acc_all_cells <- AddMetaData(object = acc_all_cells, metadata = cnsig, col.name = "copynumber")
acc_all_cells = SetIdent(object = acc_all_cells,value = "cell.type")
FeaturePlot(acc_all_cells, "copynumber",pt.size = 1, cols = c("lightblue","orange","red","darkred"),label = T,repel = T)
acc_cancer_cells <- AddMetaData(object = acc_cancer_cells, metadata = cnsig, col.name = "copynumber")
acc_cancer_cells = SetIdent(object = acc_cancer_cells,value = "cell.type")
FeaturePlot(acc_cancer_cells, "copynumber",pt.size = 1, cols = c("lightblue","orange","red","darkred"),label = T,repel = T)
CNV plot
cnv_subtypes = as.data.frame(cutree(infercnv_obj_default@tumor_subclusters[["hc"]][["HMSC"]], k = 2))
names(cnv_subtypes)[1] = "cnv.cluster"
rownames(cnv_subtypes) = rownames(cnv_subtypes) %>% gsub(pattern = "-",replacement = "\\.")
cnv_subtypes [,1] = as.character(cnv_subtypes[,1])
infercnv.observations = data.frame(fread(file = "./Data/inferCNV/infercnv.observations.txt"), row.names=1)
names_to_keep = colnames(infercnv.observations) %in% (colnames(acc1_cancer_cells) %>% gsub(pattern = "_",replacement = "\\."))
infercnv.observations = infercnv.observations[,names_to_keep]
rownames(cnv_subtypes) = rownames(cnv_subtypes) %>% gsub(pattern = "2\\.",replacement = "2_")
rownames(cnv_subtypes) = rownames(cnv_subtypes) %>% gsub(pattern = "3\\.",replacement = "3_")
acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = cnv_subtypes)
DimPlot(acc1_cancer_cells,group.by = "cnv.cluster",pt.size = 2)
original_myo_genes = c( "TP63", "TP73", "CAV1", "CDH3", "KRT5", "KRT14", "ACTA2", "TAGLN", "MYLK", "DKK3")
original_lum_genes = c("KIT", "EHF", "ELF5", "KRT7", "CLDN3", "CLDN4", "CD24", "LGALS3", "LCN2", "SLPI" )
calculate_score(dataset = all_acc_cancer_cells,myo_genes = original_myo_genes,lum_genes = original_lum_genes)
correlation of lum score and myo score: -0.45
correlation of lum score and original lum score: 1
correlation of myo score and original myo score: 1
calculate_score(dataset = acc1_cancer_cells,myo_genes = original_myo_genes,lum_genes = original_lum_genes,lum_threshold = 0,myo_threshold = 0)
Warning: Could not find TP73 in the default search locations, found in RNA assay instead
Warning: Could not find MYLK in the default search locations, found in RNA assay instead
Warning: Could not find CLDN3 in the default search locations, found in RNA assay instead
Warning: Could not find CD24 in the default search locations, found in RNA assay instead
correlation of lum score and myo score: -0.11
Only HMSC cancer cells:
HPV33_P3 = fread("./Data/HPV33_P3.txt",col.names = c("plate","reads")) %>% as.data.frame()
HPV33_P3.df = HPV33_P3 %>% mutate(
plate = gsub(x =HPV33_P3$plate, replacement = "",pattern = "_.*$")
%>% gsub(pattern = "-P",replacement = ".P")
%>% gsub(pattern = "-",replacement = "_",)
)
HPV33_P3.df = HPV33_P3.df %>% dplyr::filter(HPV33_P3.df$plate %in% colnames(acc1_cancer_cells))
rownames(HPV33_P3.df) <- HPV33_P3.df$plate
HPV33_P3.df$plate = NULL
HPV33_P2 = fread("./Data/HPV33_P2.txt",col.names = c("plate","reads")) %>% as.data.frame()
HPV33_P2.df = HPV33_P2 %>% mutate(
plate = gsub(x =HPV33_P2$plate, replacement = "",pattern = "_.*$")
%>% gsub(pattern = "plate2-",replacement = "plate2_",)
%>% gsub(pattern = "-",replacement = "\\.",)
)
HPV33_P2.df = HPV33_P2.df %>% dplyr::filter(HPV33_P2.df$plate %in% colnames(acc1_cancer_cells))
rownames(HPV33_P2.df) <- HPV33_P2.df$plate
HPV33_P2.df$plate = NULL
HPV33 = rbind(HPV33_P3.df,HPV33_P2.df)
acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = HPV33,col.name = "HPV33.reads")
FeaturePlot(acc1_cancer_cells,features = "HPV33.reads",max.cutoff = 10)
data = FetchData(object = acc1_cancer_cells,vars = "HPV33.reads")
data = data %>% mutate("0 reads" = if_else(condition = HPV33.reads == 0,true = 1,false = 0))
data = data %>% mutate("1 reads" = if_else(condition = HPV33.reads == 1,true = 1,false = 0))
data = data %>% mutate("2 reads" = if_else(condition = HPV33.reads == 2,true = 1,false = 0))
data = data %>% mutate("3-23 reads" = if_else(condition = HPV33.reads >=3 &HPV33.reads <24,true = 1,false = 0))
data = data %>% mutate("24+ reads" = if_else(condition = HPV33.reads >=24,true = 1,false = 0))
data = colSums(data[,2:ncol(data)]) %>% as.data.frame()
names(data)[1] = "count"
data = rownames_to_column(data,var = "bin")
data
ggplot(data=data, aes(x=factor(bin,levels = c("0 reads","1 reads","2 reads","3-23 reads","24+ reads")), y=count)) +
geom_bar(stat="identity", fill="steelblue") + xlab("HPV Reads")+ theme_minimal()+
geom_text(aes(label=count), vjust=-0.5, color="black", size=3.5)
hpv33_positive = HPV33 %>% dplyr::mutate(hpv33_positive = case_when(reads >= 10 ~ "positive",
reads < 10 ~ "negative")
)
hpv33_positive$reads = NULL
acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = hpv33_positive)
DimPlot(object = acc1_cancer_cells,group.by = c("hpv33_positive"),pt.size = 2)
acc1_cancer_cells$plate = acc1_cancer_cells$orig.ident
saveRDS(object = acc1_cancer_cells, file = "./Data/acc1_cancer_cells_15KnCount_V4.RDS")
library(reticulate)
acc1_cancer_cells = FindVariableFeatures(object = acc1_cancer_cells,nfeatures = nfeatures)
Calculating gene variances
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
from cnmf import cNMF
import numpy as np
nfeatures_name = r.nfeatures_name
name = 'HMSC_cNMF_'+nfeatures_name+'Kvargenes'
outdir = './Data/cNMF'
K_range = np.arange(3,10)
cnmf_obj = cNMF(output_dir=outdir, name=name)
counts_fn='./Data/cNMF/hmsc_expressionData_'+nfeatures_name+'Kvargenes.txt'
tpm_fn = counts_fn ## This is a weird case where because this dataset is not 3' end umi sequencing, we opted to use the TPM matrix as the input matrix rather than the count matrix
cnmf_obj.prepare(counts_fn=counts_fn, components=K_range, seed=14,tpm_fn=tpm_fn)
cnmf_obj.factorize(worker_i=0, total_workers=1)
cnmf_obj.combine()
cnmf_obj.k_selection_plot()
# import pickle
# f = open('./Data/cNMF/HMSC_cNMF_'+nfeatures_name+'Kvargenes/cnmf_obj.pckl', 'wb')
# pickle.dump(cnmf_obj, f)
# f.close()
from cnmf import cNMF
import pickle
nfeatures = "2K"
f = open('./Data/cNMF/HMSC_cNMF_' + nfeatures+ 'vargenes/cnmf_obj.pckl', 'rb')
cnmf_obj = pickle.load(f)
f.close()
selected_k = 4
density_threshold = 0.1
cnmf_obj.consensus(k=selected_k, density_threshold=density_threshold,show_clustering=True)
/sci/labs/yotamd/lab_share/avishai.wizel/python_envs/miniconda/envs/cnmf_env_6/lib/python3.7/site-packages/scanpy/preprocessing/_simple.py:843: UserWarning: Received a view of an AnnData. Making a copy.
view_to_actual(adata)
usage_norm, gep_scores, gep_tpm, topgenes = cnmf_obj.load_results(K=selected_k, density_threshold=density_threshold)
gep_scores = py$gep_scores
gep_tpm = py$gep_tpm
all_metagenes= py$usage_norm
plt_list = list()
for (i in 1:ncol(gep_scores)) {
top_genes = gep_scores %>% arrange(desc(gep_scores[i])) #sort by score a
top = head(rownames(top_genes),200) #take top top_genes_num
res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title =
i,silent = T,return_all = T)
plt_list[[i]] = res$plt
}
gridExtra::grid.arrange(grobs = plt_list)
canonical_pathways = msigdbr(species = "Homo sapiens", category = "C2") %>% dplyr::filter(gs_subcat != "CGP") %>% dplyr::distinct(gs_name, gene_symbol)
plt_list = list()
for (i in 1:ncol(gep_scores)) {
top_genes = gep_scores %>% arrange(desc(gep_scores[i])) #sort by score a
top = head(rownames(top_genes),200) #take top top_genes_num
res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title =
i,silent = T,return_all = T,custom_pathways = canonical_pathways)
plt_list[[i]] = res$plt
}
gridExtra::grid.arrange(grobs = plt_list)
# lum genes in metagenes
message("lum genes in metagenes")
lum genes in metagenes
for (i in 1:ncol(gep_scores)) {
top_genes = gep_scores %>% arrange(desc(gep_scores[i])) #sort by score a
top = head(rownames(top_genes),200) #take top top_genes_num
cat(paste0("metagene ",i,": "))
print(original_lum_genes[original_lum_genes %in% top])
}
metagene 1: [1] "KRT7" "SLPI"
metagene 2: [1] "CLDN4"
metagene 3: character(0)
metagene 4: [1] "KRT7" "LGALS3" "LCN2" "SLPI"
cat("\n")
# myo genes in metagenes
message("myo genes in metagenes")
myo genes in metagenes
for (i in 1:ncol(gep_scores)) {
top_genes = gep_scores %>% arrange(desc(gep_scores[i])) #sort by score a
top = head(rownames(top_genes),200) #take top top_genes_num
cat(paste0("metagene ",i,": "))
print(original_myo_genes[original_myo_genes %in% top])
}
metagene 1: [1] "KRT14" "ACTA2" "DKK3"
metagene 2: character(0)
metagene 3: character(0)
metagene 4: character(0)
cat("\n")
notch_genes = c("JAG1","JAG2","NOTCH3","NOTCH2","NOTCH1","DLL1","MYB","HES4","HEY1","HEY2","NRARP")
# notch genes in metagenes
message("myo genes in metagenes")
myo genes in metagenes
for (i in 1:ncol(gep_scores)) {
top_genes = gep_scores %>% arrange(desc(gep_scores[i])) #sort by score a
top = head(rownames(top_genes),200) #take top top_genes_num
cat(paste0("metagene ",i,": "))
print(notch_genes[notch_genes %in% top])
}
metagene 1: character(0)
metagene 2: character(0)
metagene 3: character(0)
metagene 4: character(0)
cat("\n")
plt_list = list()
for (i in 1:ncol(gep_scores)) {
top_genes = gep_scores %>% arrange(desc(gep_scores[i])) #sort by score a
top = head(rownames(top_genes),10) #take top top_genes_num
message(paste("program ",i,"top genes:"))
print(top)
}
program 1 top genes:
[1] "IGKV5-2" "NDUFB7" "LGALS1" "UBL5" "S100A6" "S100A11" "CUTA" "IFITM3" "JTB" "IL17RB"
program 2 top genes:
[1] "EGR1" "C6orf62" "JUNB" "DNAJA1" "ATF3" "IER2" "ERRFI1" "KLF6" "CDKN1A" "MTHFD2"
program 3 top genes:
[1] "RP1-128M12.3" "RP11-374F3.3" "RP11-403A21.3" "RP11-454L9.2" "AC097500.1" "RP11-463H24.4" "RP11-515O17.2"
[8] "RP11-536L3.4" "PALD1" "AL049758.2"
program 4 top genes:
[1] "LTF" "PIGR" "FMO2" "HSD17B2" "MLPH" "PRR15L" "RF00019.219"
[8] "RP11-817J15.2" "CD14" "CLU"
meta3 = FetchData(object = acc1_cancer_cells,vars = c("metagene.3"))
ggplot(meta3, aes(x=metagene.3)) +
geom_density()
meta3[,1] = as.numeric(meta3[,1])
sum(meta3[,1]>0,na.rm = T )
[1] 82
sum(meta3[,1]==0,na.rm = T )
[1] 53
larger_by = 2
message(paste("larger_by = ", larger_by))
larger_by = 2
acc1_cancer_cells = program_assignment(dataset = acc1_cancer_cells,larger_by = larger_by,program_names = colnames(all_metagenes))
selected_k =4
colors = rainbow(selected_k)
colors = c(colors,"grey")
DimPlot(acc1_cancer_cells,group.by = "program.assignment",pt.size = 2,cols =colors)
show cell cycle program:
hallmark_name = "GO_MITOTIC_CELL_CYCLE"
genesets =GSEABase::getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
var_features=acc1_cancer_cells@assays$RNA@var.features
geneIds= genesets[[hallmark_name]]@geneIds
score <- apply(acc1_cancer_cells@assays$RNA@data[intersect(geneIds,var_features),],2,mean)
acc1_cancer_cells=AddMetaData(acc1_cancer_cells,score,hallmark_name)
FeaturePlot(object = acc1_cancer_cells,features = hallmark_name)
cc_vs_program2 = FetchData(object = acc1_cancer_cells,vars = c("metagene.2",hallmark_name))
cor(cc_vs_program2[1],cc_vs_program2[2])
GO_MITOTIC_CELL_CYCLE
metagene.2 0.6127281
cnv_vs_hpv = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","hpv33_positive"))
test <- fisher.test(table(cnv_vs_hpv))
ggbarstats(
cnv_vs_hpv, cnv.cluster, hpv33_positive,
results.subtitle = FALSE,
subtitle = paste0(
"Fisher's exact test", ", p-value = ",
ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
)
)
cnv_vs_cnmf = FetchData(object = acc1_cancer_cells,vars = c("program.assignment","cnv.cluster"))
cnv_vs_cnmf = cnv_vs_cnmf %>% dplyr::filter(program.assignment == "1" |program.assignment == "2" )
test <- fisher.test(table(cnv_vs_cnmf))
ggbarstats(
cnv_vs_cnmf, program.assignment, cnv.cluster,
results.subtitle = FALSE,
subtitle = paste0(
"Fisher's exact test", ", p-value = ",
ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
)
)
cnmf_vs_hpv = FetchData(object = acc1_cancer_cells,vars = c("program.assignment","hpv33_positive"))
cnmf_vs_hpv = cnmf_vs_hpv %>% dplyr::filter(program.assignment == "1" |program.assignment == "2" )
test <- fisher.test(table(cnmf_vs_hpv))
ggbarstats(
cnmf_vs_hpv, program.assignment, hpv33_positive,
results.subtitle = FALSE,
subtitle = paste0(
"Fisher's exact test", ", p-value = ",
ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
)
)
myb_vs_cnv = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","MYB"))
myb_vs_cnv $cnv.cluster = as.character(myb_vs_cnv $cnv.cluster )
ggboxplot(myb_vs_cnv, x = "cnv.cluster", y = "MYB",
palette = "jco",
add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("1","2")))
myb_vs_hpv = FetchData(object = acc1_cancer_cells,vars = c("hpv33_positive","MYB"))
myb_vs_hpv $hpv33_positive = as.character(myb_vs_hpv $hpv33_positive )
ggboxplot(myb_vs_hpv, x = "hpv33_positive", y = "MYB",
palette = "jco",
add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("positive","negative")))+ stat_summary(fun.data = function(x) data.frame(y=15, label = paste("Mean=",round(mean(x),digits = 2))), geom="text") +ylab("log2(MYB)")
plate_1_new = subset(acc1_cancer_cells,subset = orig.ident == "ACC.plate2")
myb_vs_hpv = FetchData(object = plate_1_new,vars = c("hpv33_positive","MYB"))
myb_vs_hpv $hpv33_positive = as.character(myb_vs_hpv $hpv33_positive )
# myb_vs_hpv$MYB = 2**myb_vs_hpv$MYB
p = ggboxplot(myb_vs_hpv, x = "hpv33_positive", y = "MYB",
palette = "jco",
add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("positive","negative")))+ stat_summary(fun.data = function(x) data.frame(y=15, label = paste("Mean=",round(mean(x),digits = 2))), geom="text") +ylab("log2(MYB)")+ggtitle("ACC.plate2")
plate_2_new = subset(acc1_cancer_cells,subset = orig.ident == "ACC1.P3")
myb_vs_hpv = FetchData(object = plate_2_new,vars = c("hpv33_positive","MYB"))
myb_vs_hpv $hpv33_positive = as.character(myb_vs_hpv $hpv33_positive )
# myb_vs_hpv$MYB = 2**myb_vs_hpv$MYB
p+
ggboxplot(myb_vs_hpv, x = "hpv33_positive", y = "MYB",
palette = "jco",
add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("positive","negative")))+ stat_summary(fun.data = function(x) data.frame(y=15, label = paste("Mean=",round(mean(x),digits = 2))), geom="text") +ylab("log2(MYB)")+ggtitle("ACC1.P3")
Warning in wilcox.test.default(c(4.76973518359805, 7.43610462068823, 4.9899976998433, :
cannot compute exact p-value with ties
hpvReads_vs_myb = FetchData(object = acc1_cancer_cells,vars = c("HPV33.reads","MYB"))
corr = cor(hpvReads_vs_myb$HPV33.reads,hpvReads_vs_myb$MYB)
print("correlation of MYB abd hpv33_reads:" %>% paste(corr %>% round(digits = 2)))
acc1_cancer_cells_data = acc1_cancer_cells@assays[["RNA"]]@data %>% as.data.frame()
acc1_cancer_cells_data = cor(acc1_cancer_cells_data)
annotation = FetchData(object = acc1_cancer_cells,vars = c("orig.ident"))
colors <- c(seq(-1,1,by=0.01))
my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
(n = length(colors)-3), "red")
pht1 = pheatmap(acc1_cor,annotation_col = annotation,fontsize = 6,breaks = colors, color = my_palette,show_colnames = F,show_rownames = F)
acc1_cancer_cells_data = acc1_cancer_cells@assays[["RNA"]]@data %>% as.data.frame()
acc1_cancer_cells_data = cor(acc1_cancer_cells_data)
annotation = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster"))
colors <- c(seq(-1,1,by=0.01))
my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
(n = length(colors)-3), "red")
pht1 = pheatmap(acc1_cor,annotation_col = annotation,fontsize = 6,breaks = colors, color = my_palette,show_colnames = F,show_rownames = F)
cnv_vs_plate = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","orig.ident"))
test <- fisher.test(table(cnv_vs_plate))
ggbarstats(
cnv_vs_plate, cnv.cluster, orig.ident,
results.subtitle = FALSE,
subtitle = paste0(
"Fisher's exact test", ", p-value = ",
ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
)
)
# creat colours for each group
cnv_vs_plate = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","orig.ident","genes_detected"))
cnv_vs_plate$genes_detected[cnv_vs_plate$genes_detected>8000] = 8000
rownames(cnv_vs_plate) = rownames(cnv_vs_plate) %>% gsub(pattern = "_",replacement = "\\.")
cnv_vs_plate$cnv.cluster = as.character(cnv_vs_plate$cnv.cluster)
cnv_vs_plate = cnv_vs_plate %>% dplyr::rename(plate = "orig.ident")
cnv_vs_plate = cnv_vs_plate %>% arrange(plate,genes_detected)
infercnv.observations2_ordered <- infercnv.observations2[match(rownames(cnv_vs_plate), rownames(infercnv.observations2) ),]
annoCol <- list(plate = c(ACC1.P3 = "red",ACC.plate2 = "green"),cnv.cluster =c("1"="green","2" = "red"))
pheatmap(infercnv.observations2_ordered,cluster_cols = F,cluster_rows = F, show_rownames = F,show_colnames = F, breaks = breaks,color = colorRampPalette(rev(c("darkred", "white", "darkblue")))(15),annotation_row = cnv_vs_plate,annotation_colors = annoCol)
# creat colours for each group
cnv_vs_plate = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","orig.ident"))
rownames(cnv_vs_plate) = rownames(cnv_vs_plate) %>% gsub(pattern = "_",replacement = "\\.")
cnv_vs_plate$cnv.cluster = as.character(cnv_vs_plate$cnv.cluster)
cnv_vs_plate = cnv_vs_plate %>% dplyr::rename(plate = "orig.ident")
annoCol <- list(plate = c(ACC1.P3 = "red",ACC.plate2 = "green"),cnv.cluster =c("1"="green","2" = "red"))
pheatmap(infercnv.observations2,cluster_cols = F,cluster_rows = F, show_rownames = F,show_colnames = F, breaks = breaks,color = colorRampPalette(rev(c("darkred", "white", "darkblue")))(15),annotation_row = cnv_vs_plate,annotation_colors = annoCol)
data = FetchData(object = acc1_cancer_cells,vars = c("metagene.1","metagene.2"))
# data$genes_detected = min_max_normalize(data$genes_detected)
annotation = FetchData(object = acc1_cancer_cells,vars = c("genes_detected","program.assignment","orig.ident"))
annotation = annotation %>% dplyr::rename(plate = "orig.ident")
annotation$genes_detected[annotation$genes_detected>8000] = 8000
annotation = annotation %>% arrange(plate,genes_detected)
data <- data[match(rownames(annotation), rownames(data) ),]
annoCol <- list(program.assignment = c(metagene.1 = "red",metagene.2 = "grey96",metagene.3 = "yellow",metagene.4 = "blue","NA" = "grey"),plate = c(ACC1.P3 = "red",ACC.plate2 = "green"))
colors <- c(seq(-1,1,by=0.01))
my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
(n = length(colors)-3), "red")
pheatmap(data,cluster_cols = F,cluster_rows = F,show_rownames = F,color = my_palette, breaks = colors,annotation_row = annotation,annotation_colors = annoCol)
Warning in grSoftVersion() :
unable to load shared object '/usr/local/lib/R/modules//R_X11.so':
libXt.so.6: cannot open shared object file: No such file or directory
plate_3 = cnmf_vs_plate %>% dplyr::filter((program.assignment == 1 & orig.ident == "ACC1.P3") ) %>% rownames()
plate_2 = cnmf_vs_plate %>% dplyr::filter((program.assignment == 2 & orig.ident == "ACC.plate2")) %>% rownames()
cells = list(ACC1.P3 = plate_3,ACC.plate2 = plate_2)
exceptions_plt = DimPlot(object = acc1_cancer_cells, cells.highlight = cells, cols.highlight = c("green","red"), cols = "gray", order = TRUE,pt.size = 2,sizes.highlight = 2,combine = F)
Warning messages: 1: No Python documentation found for ‘cnmf_obj.prepare’. 2: No Python documentation found for ‘cnmf_obj.prepare’.
colors = c("green","red","cyan","purple","grey")
program.assignment_plt = DimPlot(acc1_cancer_cells,group.by = "program.assignment",pt.size = 2,cols = colors,combine = F)
metagene.1_plt <- FeaturePlot(object = acc1_cancer_cells,features = c("metagene.1"),combine = F)
metagene.2_plt = FeaturePlot(object = acc1_cancer_cells,features = c("metagene.2"),combine = F)
lst = list(exceptions = exceptions_plt, program.assignment = program.assignment_plt,metagene.1 = metagene.1_plt,metagene.2 = metagene.2_plt)
for (i in 1:length(lst)) {
cat("### ",names(lst)[i]," \n")
print(
lst[[i]]
)
plot.new()
dev.off()
cat(' \n\n')
}
[[1]] Warning in grSoftVersion() : unable to load shared object ‘/usr/local/lib/R/modules//R_X11.so’: libXt.so.6: cannot open shared object file: No such file or directory
[[1]]
[[1]]
[[1]]
NA